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ABSTRACT 

The  vector  held  relations  and  mobile  charge  thermo¬ 
dynamics  that  form  a  complete  and  self-consistent  classical 
model  of  the  electromagnetic  interaction  are  presented.  The 
charge  thermodynamics  includes  an  original  treatment  of 
heat  how  between  ideal  Fermi  gases  that  is  derived  from 
their  heat  capacities.  An  original  discretization  scheme 
based  on  the  properties  of  Delaunay  and  Voronoi  meshes 
is  also  presented.  This  new  scheme  allows  the  held  equa¬ 
tions  to  be  solved  self-consistently  with  the  highly  nonlinear 
charge  transport  equations,  producing  the  fully  coupled  dy¬ 
namics  of  full  wave  vector  helds,  mobile  charge  densities, 
as  well  as  mobile  charge  and  crystal  lattice  temperatures. 
Linear  and  nonlinear  lossy  transmission  lines  are  used  to 
demonstrate  the  simulator. 

1.  INTRODUCTION 

Decades  of  research  in  computational  electromagnet¬ 
ics  and  electronics  have  led  to  both  sophisticated  full  wave 
held  solvers  and  advanced  semiconductor  device  simulators. 
However,  both  types  of  computation  typically  incorporate 
certain  simplifying  assumptions  that  limit  their  applica¬ 
tions.  For  example,  full  wave  held  solvers  commonly  ap¬ 
proximate  charge  transport  as  a  linear  process  in  which  the 
mobile  charge  hux  is  directly  proportional  to  the  electric 
held.  On  the  other  hand,  while  active  semiconductor  de¬ 
vice  simulators  offer  a  much  more  detailed  description  of 
nonlinear  charge  dynamics,  they  typically  treat  the  helds 
with  the  electrostatic  approximation.  Although  there  are 
many  important  applications  for  these  two  distinctly  differ¬ 
ent  treatments  of  the  electromagnetic  interaction,  current 
technological  trends  have  created  a  need  to  unify  the  two 
approaches. 

Advances  in  high  performance  microelectronics  have 
produced  smaller  component  sizes,  higher  circuit  densities, 
and  increased  operating  frequencies.  These  advances  have 
led  to  full  wave  phenomena  that  affect  the  behavior  of  ac¬ 
tive  components  and  circuits.  For  example,  at  clock  speeds 
above  a  few  hundred  megahertz,  surface  waves,  propagation 
delays,  as  well  as  unintended  radiation,  causing  crosstalk 
between  densely  integrated  components,  can  become  im¬ 
portant.  At  still  higher  frequencies,  e.g.  in  the  millimeter 
wave  range,  there  may  even  be  a  qualitative  shift  in  the  way 
individual  active  components  behave.  As  fields  fluctuate  at 
rates  comparable  to  mobile  charge  momentum  and  energy 
relaxation  rates,  an  active  device  can  behave  less  as  a  non- 
linearly  component  that  manipulates  a  signal  and  more  like 
a  linear  passive  element  that  stores  electromagnetic  energy. 


Aside  from  unintended  consequences,  other  advanced 
technologies  exploit  the  tight  nonlinear  coupling  between 
electromagnetic  radiation  and  mobile  charges.  For  example, 
a  possible  ultrahigh  frequency  modulation  scheme  for  laser 
diodes  is  based  on  heating  degenerate  two-dimensional  gases 
of  mobile  charges  in  a  quantum  well  with  terahertz  radia¬ 
tion  (Li  and  Ning  2000).  Another  new  technology  uses  hot 
electron  induced  real  space  transfer  in  the  channel  of  a  short 
gate  length  HEMT  to  emit  terahertz  radiation,  providing  a 
useful  compact  source  for  radiation  in  this  frequency  range 
(Knap  et  al.  2004).  Simulating  nonlinear  charge  transport 
with  full  wave  electromagnetics  is  required  to  advance  these 
new  technologies  as  well  as  to  understand  the  full  wave  par- 
asitics  in  high  speed  microelectronics. 

The  need  for  full  wave  treatment  of  active  devices  has 
led  to  work  on  combining  full  wave  field  solvers  and  semicon¬ 
ductor  device  simulators.  While  virtually  all  these  efforts 
treat  the  fields  using  the  Finite  Difference  Time  Domain 
(FDTD)  method  developed  by  Yee  (Yee  1966),  they  differ 
in  the  levels  of  detail  with  which  they  treat  active  devices. 
Some  approaches  replace  the  active  devices  with  lumped 
parameter  equivalent  circuits  which  impose  field  controlled 
current  sources  on  certain  edges  of  the  FDTD  mesh  (Sui 
et  al.  1992),  (Piket-May  et  al.  1994),  (Ciampolini  et  al. 
1996).  Other  approaches  treat  carrier  transport  directly  by 
coupling  the  FDTD  field  equations  to  the  Boltzmann  equa¬ 
tion,  solved  classically  with  drift-diffusion  or  hydrodynamic 
theory  (Sano  and  Shibata  1990),  (Alsunaidi  et  al.  1996), 
(Mavahhedi  and  Abdipour  2006)  or  semi-classically  with  the 
Monte  Carlo  method  (El-Ghazaly  et  al.  1990),  (Goodnick 
et  al.  1995). 

Although  different  full  wave  simulations  of  active  de¬ 
vices  treat  charge  transport  differently,  their  common  use 
of  FDTD  field  discretization  gives  them  certain  similari¬ 
ties.  For  example,  FDTD  requires  a  structured  mesh  and 
uses  explicit  time  stepping.  Explicit  time  stepping  can  be 
performed  without  factoring  a  large  Jacobian  matrix  to  de¬ 
termine  the  rotational  fields.  On  the  surface,  this  is  a  great 
advantage,  allowing  the  solution  of  problems  with  tens  of 
millions  of  unknowns  (Piket-May  et  al.  1994).  However,  this 
advantage  is  undone,  at  least  in  part,  by  the  highly  non¬ 
linear  coupling  of  the  fields  and  the  mobile  charge  fluxes. 
This  tight  coupling  imposes  severe  discretization  restric¬ 
tions  that  can  require  dense  meshes,  refined  in  terms  of 
the  Debye  length,  to  represent  the  Boltzmann  equation  (Al¬ 
sunaidi  et  al.  1996),  (Mavahhedi  and  Abdipour  2006)  and 
very  short  time  steps  determined  by  the  Courant  condition. 
With  these  discretization  requirements,  solving  a  relatively 
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small  problem  domain  can  require  millions  of  solution  vari¬ 
ables  solved  repeatedly  for  tens  of  thousands  of  time  steps. 
Finally,  the  “leapfrog”  scheme  FDTD  uses  to  step  through 
time  means  that  at  no  moment  during  the  simulation  are 
the  electric  and  magnetic  fields  and  the  mobile  charge  fluxes 
self-consistent.  In  addition  to  adding  to  stability  problems, 
lack  of  self-consistency  means  that  the  highly  nonlinear  ex¬ 
change  of  energy  between  the  electromagnetic  fields  and  the 
mobile  charges  is  never  fully  represented. 

This  paper  presents  a  self-consistent  classical  model 
of  the  electromagnetic  interaction  suitable  for  full  wave 
simulation  of  active  semiconductor  devices.  Electromag¬ 
netic  fields  are  discretized  using  a  new  scheme  called  De- 
launay/Voronoi  Surface  Integration  (DVSI)  which  allows 
unstructured  tetrahedral  meshes.  The  mobile  charge  and 
energy  fluxes  are  derived  from  moments  of  the  Boltzmann 
equation  and  incorporate  a  new  treatment  of  heat  flow  de¬ 
rived  from  the  heat  capacity  of  Fermi  gases.  The  field  equa¬ 
tions  are  solved  self-consistently  with  the  charge  and  en¬ 
ergy  transport  equations.  This  greatly  improves  the  sta¬ 
bility  of  the  simulator  so  that  discretization  is  determined 
solely  by  variations  in  the  solution  variables.  That  is,  the 
mesh  need  only  be  refined  enough  to  adequately  represent 
field  wavelengths  and  rapid  changes  in  charge  densities,  e.g. 
dopant  junctions,  heterojunctions,  Schottky  barriers.  Like¬ 
wise,  time  steps  need  only  by  small  enough  to  represent 
dynamic  variations  in  the  solution  variables.  This  permits 
adaptive  time  stepping  where  time  steps  may  vary  by  several 
orders  of  magnitude  during  the  course  of  a  simulation.  To 
demonstrate  the  model,  certain  linear  and  nonlinear  lossy 
transmission  line  simulations  are  presented. 

2.  FULL  WAVE  NONLINEAR 
CHARGE  TRANSPORT 

Maxwell’s  vector  field  theory  and  his  kinetic  theory  of 
gases  form  a  complete  and  self-consistent  model  of  the  elec¬ 
tromagnetic  interaction.  The  model  consists  of  Maxwell’s 
static  and  dynamic  field  equations  as  well  as  the  equations 
for  charge  and  energy  conservation.  Charge  and  energy  con¬ 
servation  can  be  derived  from  moments  of  the  Boltzmann 
equation,  and  the  fluxes  that  appear  in  these  conservation 
laws  can  also  be  obtained  from  moments  of  the  Boltzmann 
equation  with  special  attention  payed  to  the  thermodynam¬ 
ics  of  ideal  gases. 

2.1  Maxwell’s  Field  Equations 

The  behavior  of  electromagnetic  fields  is  completely  de¬ 
termined  by  Maxwell’s  equations. 


V  •  eE  =  p 

(1) 

X7  -pH  =  0 

(2) 

deE 

(3) 

V  x  H  =  J  -\ — — — 
at 

<^> 

cc 

1 

II 

X 

> 

(4) 

where  (1)  and  (2)  are  Gauss’s  electrostatic  and  magneto¬ 
static  laws  respectively,  (3)  is  Ampere’s  law,  (4)  is  Faraday’s 


law,  e  is  the  dielectric  constant,  p  is  the  magnetic  perme¬ 
ability,  E  is  the  electric  field,  H  is  the  magnetic  field,  p  is 
the  charge  density,  and  J  is  the  mobile  charge  flux.  The 
charge  density  is  determined  locally  by  the  concentration 
of  mobile  electrons  n,  mobile  holes  p,  ionized  donor  impu¬ 
rities  N+,  and  ionized  acceptors  1VJ.  The  charge  flux  is 
determined  locally  by  mobile  electron  and  hole  fluxes.  The 
charge  densities  and  fluxes  are  obtained  from  the  transport 
model. 

To  find  the  electric  and  magnetic  fields  that  satisfy  the 
four  Maxwell’s  equations,  the  electric  field  can  be  expressed 
in  terms  of  vector  and  scalar  potentials,  A  and  if)  respec¬ 
tively. 

/  d  A  \ 

E  =  —  (  +  V'i/’J  =  Erot  +  Ejrr  (5) 

Since  the  curl  of  a  gradient  is  always  zero,  Vi/;  represents 
a  conservative,  irrotational  electric  field.  To  make  the  two 
electric  field  components  functionally  orthogonal,  a  zero  di¬ 
vergence  is  assigned  to  the  divergence  of  the  vector  poten¬ 
tial’s  displacement. 


This  is  equivalent  to  the  Coulomb  gauge  for  a  homogeneous 
medium.  It  allows  the  time  derivative  of  the  vector  potential 
to  be  viewed  as  an  entirely  rotational  electric  field. 

From  Faraday’s  law,  the  magnetic  flux  density  is  given 
by  the  curl  of  the  vector  potential. 

pH  =  V  x  A 

Since  the  divergence  of  a  curl  is  always  zero,  Gauss’s  mag¬ 
netostatic  law  is  automatically  satisfied.  This  means  that 
the  magnetic  field  is  treated  as  entirely  rotational,  which  is 
valid  for  homogeneous  permeability.  For  nonhomogeneous 
p,  e.g.  ferrites,  the  magnetic  field  should  be  expressed  in 
terms  of  rotational  and  conservative  components,  similar  to 
(5).  However,  this  paper  only  considers  homogeneous  per¬ 
meability. 

2.2  Classical  Charge  Transport 

For  classical  systems,  charges  can  be  treated  as  particles 
with  their  behavior  determined  by  the  Boltzmann  equation. 
A  relatively  efficient  way  to  solve  this  equation  treats  the 
charges  as  ideal  Fermi  gases.  An  ideal  Fermi  gas  in  three 
dimensions  is  spherically  symmetric  in  momentum  space, 
and  its  distribution  in  energy  is  determined  by  its  chemical 
potential  and  temperature  according  to  Fermi-Dirac  statis¬ 
tics. 


where  fn  is  the  occupation  probability  of  an  electron  in 
the  conduction  band  of  a  semiconductor,  E  is  the  electron 
energy,  Fn  is  its  chemical  potential  (quasi-Fermi  level),  and 
Tn  is  the  electron  temperature.  An  analogous  expression  ex¬ 
ists  for  the  distribution  function  of  positively  charged  holes 
in  the  valence  band.  The  effective  mass  approximation  for 
band  structure  can  be  used  to  integrate  /„  over  all  available 
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momentum  states  and  define  the  local  electron  density. 


V2ml/2  f°°  VE-EC 

n  =  - ~ —  /  - 7 - r - ah 

TT2h‘  J ec  exp  ^  Ekt"  )  +  1 

=  Nc(kT)3/2F1/2(Vn ) 

where  mn  is  the  effective  mass  for  electrons,  Ec  is  the  con¬ 
duction  band  edge,  F,,  is  the  Fermi  integral  of  degree  i,  and 
Vn  =  (En  —  Ec)/{kTn).  A  similar  expression  exists  for  the 
density  of  holes  in  the  valence  band.  Fermi-Dirac  statistics 
is  also  used  to  calculate  the  density  of  ionized  dopant  impu¬ 
rities.  The  electron  energy  density  is  obtained  in  a  similar 
manner. 

En  =  Nc(kT)5/2F3/2(Vn) 

The  classical  charge  transport  model  is  obtained  from 
moments  of  the  Boltzmann  equation,  expressed  in  terms  of 
the  relaxation  time  and  effective  mass  approximations  (Hess 
2000).  For  electrons,  the  moments  have  the  following  form. 

hJo{fn  +  T^h  =  T-F.Vkfn-rv.Vfn]idk  (6) 

where  r  is  the  momentum  relaxation  time,  v  is  the  electron 
velocity,  and  k  is  its  momentum  vector.  The  total  force 
F  acting  on  the  electron  includes  both  the  Coulomb  and 
Lorentz  forces. 


F  =  —q(E  +  v  X  B) 


Different  moments  are  obtained  with  different  choices  of  the 
operator  O  acting  on  the  Boltzmann  equation. 

Electron  conservation  is  obtained  by  setting  0  =  1  and 
integrating  over  all  fc-states  to  obtain  the  zeroth  moment. 

dn 

— "se/  =  V  •  J  n  +  Un 

at 

where  Jn  is  the  electron  flux  and  U„  is  the  net  recombi¬ 
nation  rate.  This  paper  considers  only  the  Shockley-Hall- 
Read  recombination  mechanism  (Hess  2000).  The  electron 
flux  J n  is  obtained  by  setting  O  =  v  and  integrating  over 
all  k  states  to  yield  the  drift-diffusion  flux. 


Jn  +  T 


dJn 

dt 
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where  f  is  an  average  momentum  relaxation  lifetime,  jin  is 
the  electron  mobility,  B  =  fiH  is  the  magnetic  flux  density, 
and  M  is  a  tensor  that  generates  the  vortex  flux  compo¬ 
nents  produced  by  the  Lorentz  force. 
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The  drift-diffusion  flux  (7)  is  discretized  using  the 
Scharfetter-Gummel  method  (Scharfetter  and  Gummel 
1969).  The  resulting  net  flux  along  a  mesh  edge  of  length 
L  connecting  vertices  1  and  2  is  given  by  the  following. 


+  f 


QJ 1=2 
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dt 


2  F3/2  kT 

3  F1/2  q 

HnM  Nc 

i  +  nlB2^x 
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Conservation  of  electron  energy  can  also  be  derived 
from  (6).  The  energy  continuity  equation  is  obtained  by 
setting  O  =  E  and  integrating  over  all  fc-states  to  obtain 
the  second  moment. 

=  qE .  Jn  +  v  •  ^ot  +  uE 

at 

where  S ^ot  is  the  total  electron  energy  flux  and  Ue  is  the  net 
loss  of  energy  to  both  recombination  and  phonon  scattering. 

UE  =  (E^)Un  +  (  kT-~kT^ 

U 1/2  \  Tn 

where  {E^ln}  is  the  average  kinetic  energy,  Tiat  is  the  tem¬ 
perature  of  the  semiconductor  crystal,  and  Tn  is  the  energy 
relaxation  lifetime. 

The  components  of  total  energy  flux  are  known  from 
the  thermodynamics  of  ideal  Fermi  gases.  In  drift-diffusion, 
when  a  charge  moves  through  space,  it  is  removed  from  one 
Fermi  gas  with  a  particular  chemical  potential  and  tem¬ 
perature,  and  it  is  transferred  to  another  Fermi  gas  with 
a  different  chemical  potential  and  temperature.  Obviously, 
the  internal  (kinetic)  energy  of  the  transferred  charges  is 
part  of  the  total  energy  exchanged  between  the  gases.  If 
the  gases  have  different  chemical  potentials  with  respect  to 
their  ground  states,  then  there  is  also  a  certain  amount  of 
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chemical  work  that  also  counts  as  an  energy  transfer.  Fi¬ 
nally,  if  the  gases  have  different  temperatures,  then  trans¬ 
ferred  charges  must  be  heated  or  cooled  to  the  temperature 
of  the  receiving  gas.  This  thermalization  process  has  the 
form  of  heat  flow  between  the  two  Fermi  gases. 

Setting  O  =  vE  in  (6)  and  integrating  over  all  fc-states 
is  the  third  moment  and  accounts  for  both  the  internal  en¬ 
ergy  and  chemical  work  components  of  the  total  flux.  Ap¬ 
plying  the  Scharfetter-Gummel  method  then  results  in  the 
following  discrete  expression  for  these  components  of  energy 
flux  between  neighboring  gases  1  and  2. 
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Transport  simulators  based  on  both  drift-diffusion  and 
hydrodynamics  treat  heat  flow  between  gases  according  to 
the  Stratton  model  (Stratton  1962),  which  expresses  it  in 
terms  of  a  thermal  conductivity  of  the  charges  and  a  tem¬ 
perature  gradient.  However,  the  heat  required  to  thermalize 
charges  moving  from  one  gas  to  another  can  be  more  pre¬ 
cisely  defined  in  terms  of  the  heat  capacity  of  Fermi  gases. 


Cv  =  (E  —  Fn) 


d 

d(kTn ) 


The  net  heat  flow  can  then  be  obtained  by  defining  the 
operator  in  (6)  as  O  =  vCy,  evaluating  the  temperature 
derivatives  of  the  distribution  function  /„,  integrating  over 
all  fc-states,  applying  the  Scharfetter-Gummel  discretization 
scheme,  and  integrating  the  exchanged  charges  from  their 
initial  to  final  temperatures.  The  resulting  discretized  net 
heat  flow  is  then  given  by  the  following. 

001^2 

cl— >2  1  -  heat, net  _  rnl^2  _  cl^2 
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The  vector  field  equations  in  Section  2.1  along  with 
the  charge  and  energy  fluxes  in  Section  2.2  form  a  physi¬ 
cally  self-consistent  classical  model  of  the  electromagnetic 
interaction.  When  solved  simultaneously,  they  show  how 
dynamic  fields  and  mobile  charges  interact  within  an  ideal¬ 
ized  semiconductor  that  acts  as  a  perfect  heat  sink  always 
remaining  at  room  temperature.  However,  real  semicon¬ 
ductors  have  finite  thermal  conductivities  and  heat  capaci¬ 
ties  that  result  in  lattice  heating  as  current  passes  through 
them.  To  account  for  this,  an  additional  equation  balanc¬ 
ing  the  production  of  lattice  energy  by  Joule’s  heat  with  its 
dissipation  by  thermal  diffusion  is  also  solved. 

-psCp^j^-  =  V  •  KVTiat  -  -E applied  •  q{J p  ~  J  n) 

where  ps  is  the  semiconductor’s  density,  Cp  is  its  specific 
heat,  and  k  is  its  thermal  conductivity. 

2.4  Time  Stepping 

For  general  input  excitations,  including  sharp  large  sig¬ 
nal  pulses,  mobile  charge  densities  and  temperatures  can 
initially  fluctuate  very  rapidly  (femtosecond  time  scales). 
Later  these  fluctuations  dampen  as  the  system  approaches 
a  new  steady  state.  The  damping  process  can  take  a  rela¬ 
tively  long  time,  reaching  steady  state  after  tens  or  hundreds 
of  nanoseconds.  Spanning  this  wide  dynamic  range  requires 
adaptive  time  stepping. 

All  dynamic  quantities  in  this  model  are  treated  as 
piece- wise  linear  in  time,  and  their  time  derivatives  are 
evaluated  with  the  fully  implicit  backward  Euler  method. 
Therefore,  any  solution  variable  0  at  time  step  i  is  approx¬ 
imated  by  the  first  two  terms  of  a  generally  infinite  series. 

0i  =  0i- 1  +  0*-lAt  +  -0i_iAt2  +  ■  ■  ■ 

The  third  term  in  the  series  serves  as  a  measure  of  0’s  devia¬ 
tion  from  linearity  and  as  a  useful  criterion  for  determining 
the  next  time  step.  When  0,;  is  computed,  simple  estimates 
of  its  time  derivatives  are  calculated. 
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After  a  maximum  relative  error  emax  is  chosen,  the  second 
derivative  can  then  be  used  to  determine  the  next  appro¬ 
priate  time  step. 
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where  S^J  are  the  same  energy  fluxes  defined  in  (9)  and 
SlZ3Tn  j  are  th°se  same  fluxes  calculated  with  both  gases 
held  at  Tnj.  Likewise,  J are  the  same  charge  fluxes 
defined  in  (8)  and  are  those  same  fluxes  with  both 

gases  at  Tnj.  Analogous  discrete  expressions  exist  for  the 
net  mobile  hole  charge  and  energy  fluxes,  Jp  and  S'*0*  re¬ 
spectively. 

2.3  Lattice  Heating 


Using  emax  =  0.001,  following  this  procedure  for  the  elec¬ 
tron  and  hole  chemical  potentials  and  temperatures  as  well 
as  for  the  lattice  temperature,  and  using  the  smallest  com¬ 
puted  Af,_|_i  as  the  next  time  step  has  been  an  effective  way 
to  capture  the  wide  dynamic  range  that  electron  devices  can 
exhibit. 

3.  DELAUNAY /VORONOI 
SURFACE  INTEGRATION 
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To  express  the  equations  in  Section  2  so  that  they  can 
be  solved  simultaneously,  the  useful  properties  of  a  Delau¬ 
nay  mesh  (Delaunay  1934)  representing  the  problem  domain 
can  be  employed.  A  contiguous  set  of  tetrahedrons  satisfies 
the  Delaunay  criterion  if  none  of  the  vertices  lie  within  any 
of  the  tetrahedrons’  circumspheres.  When  this  is  true,  a  re¬ 
ciprocal  mesh,  whose  vertices  are  the  circumcenters  of  the 
tetrahedrons,  can  be  computed.  The  reciprocal  mesh  is  the 
dual  of  the  primary  Delaunay  mesh,  and  it  consists  of  con¬ 
tiguous  Voronoi  polyhedrons  that  fill  the  problem  domain. 
Figure  1  shows  a  set  of  Delaunay  tetrahedrons  that  share  a 
common  primary  edge  and  the  Voronoi  polygon  that  corre¬ 
sponds  to  that  edge.  The  edge  is,  by  definition,  perpendic- 


Fig.  1:  Delaunay  tetrahedrons  sharing  a  common  edge  and 
the  Voronoi  polygon  facet  formed  by  their  circumcenters. 

ular  to  the  polygon. 

The  relationship  between  the  primary  edge  and  the  dual 
polygon  provides  an  easy  way  to  represent  Ampere’s  law. 
A  projection  of  an  electric  field  is  assigned  to  each  primary 
edge,  and  projections  of  magnetic  fields  are  assigned  to  the 
dual  edges  that  form  the  perimeter  of  the  Voronoi  polygon. 
The  value  of  the  electric  field  Ei  on  a  primary  edge  is  de¬ 
fined  by  integrating  (3)  over  the  area  of  the  polygon  Ai  and 
applying  the  Stokes  curl  theorem. 

/  (J + °w)  da  =  J  V  x  H  '  da  =  j>  H  '  dl 

(ji  +  Ai  —  HjLj  (10) 

In  this  way,  the  electric  field  projection  on  each  primary 
edge  is  defined  in  terms  of  magnetic  field  projections  on 
dual  edges. 

The  relationship  between  Delaunay  meshes  and  their 
Voronoi  reciprocals  also  offers  a  convenient  way  to  define 


the  magnetic  field  projections  assigned  to  dual  edges.  Fig¬ 
ure  2  shows  that  every  dual  edge  is  perpendicular  to  a  pri¬ 
mary  triangle.  The  projection  of  the  magnetic  field  Hj  on 


Fig.  2:  Dual  edges  are  perpendicular  to  primary  Delaunay 
triangles. 

a  dual  edge  is  defined  by  integrating  (4)  over  the  area  of  its 
corresponding  primary  triangle  Aj  and  applying  the  Stokes 
curl  theorem. 

—  [  —  ■  da  =  [  V  x  E  ■  da  =  (f  E  ■  dl 

J  Aj  Ot  JA.  JdA. 

=  ai) 

i 

Applying  this  procedure  to  each  dual  edge  and  combining 
the  results  with  the  equations  for  the  electric  field  projec¬ 
tions  (10)  produces  a  set  of  linearly  independent  equations 
that  can  be  solved  to  determine  the  electric  and  magnetic 
fields  throughout  the  problem  domain. 

In  addition  to  providing  convenient  ways  to  represent 
the  curl  operators  in  Ampere’s  and  Faraday’s  laws,  the  De¬ 
launay  and  Vornoi  meshes  can  be  used  to  represent  the  di¬ 
vergence  operators  that  appear  in  electrostatics  as  well  as 
charge  and  energy  conservation.  Referring  back  to  Figure  1, 
each  primary  edge  can  be  associated  with  a  Voronoi  poly¬ 
gon.  As  shown  in  Figure  3,  each  primary  vertex  belongs  to 
a  number  of  primary  edges,  all  of  which  have  Voronoi  poly¬ 
gons  that  fit  together  contiguously  to  form  a  polyhedron 
enclosing  the  vertex.  To  evaluate  the  divergence  of  a  flux, 
such  as  the  electric  displacement  in  Gauss’s  electrostatic 
law,  it  is  first  integrated  over  the  volume  of  the  polyhedron. 
By  the  divergence  theorem,  the  volume  integral  of  the  diver¬ 
gence  is  equivalent  to  the  integral  of  the  flux  over  the  surface 
area  of  the  polyhedron.  This  surface  integral  is  evaluated 
by  simply  visiting  each  primary  edge  to  which  the  vertex 
belongs,  evaluating  the  flux  on  the  edge,  and  multiplying  it 
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Fig.  3:  The  Voronoi  polygons  for  all  the  primary  edges  con¬ 
taining  a  particular  vertex  form  a  faceted  polyhedron  that 
encloses  the  vertex.  This  polyhedron  is  used  in  the  box  in¬ 
tegration  method  to  represent  divergences  away  from  the 
vertex. 


by  the  area  of  its  associated  dual  polygon. 

/  V  •  eEdv  =  /  eE-da  =  y' 

JVi  JdVi  j  Lij 

The  procedure  is  called  the  box  integration  method  and  is 
widely  used  in  semiconductor  device  simulation. 

To  summarize,  DVSI  is  used  to  represent  Ampere’s  and 
Faraday’s  laws  on  the  Delaunay  mesh,  and  box  integration 
is  used  to  represent  Gauss’s  electrostatic  law,  the  continu¬ 
ity  of  electrons  and  holes,  the  conservation  of  electron  and 
hole  energies,  as  well  as  the  conservation  of  lattice  energy. 
The  result  is  a  set  of  highly  nonlinear  coupled  differential 
equations  which  is  solved  simultaneously  with  the  full  New¬ 
ton  method.  The  solutions  show  the  full  wave  behavior  of 
the  electromagnetic  fields  coupled  to  the  nonlinear  thermo¬ 
dynamics  of  the  mobile  charges  at  each  moment  of  time 
considered  during  the  simulation. 


is  treated  as  uniformly  doped  (1015  cm-3  p-type)  silicon. 
The  line  is  excited  with  Gaussian  time  domain  input  sig¬ 
nals,  and  the  £>21  scattering  parameters  are  extracted  from 
Fourier  analyses  of  the  simulated  input  and  output  voltages 
and  currents.  This  simple,  highly  linear  test  case  permits 
accurate  analytic  approximations  for  the  S-parameters  de¬ 
rived  from  the  telegraphist’s  equations  (Ramo  et  al.  1984). 
Comparing  the  analytic  S-parameters  with  the  simulated 
values  is  a  useful  test  of  both  the  charge  transport  model 
and  the  field  solutions. 

Figure  5  compares  the  derived  S-parameters  with  their 
simulated  values,  computed  with  both  full  wave  electromag¬ 
netics  and  with  the  quasi-static  approximation.  The  full 


Fig.  5:  Simulated  S-parameters  for  the  lossy  transmission 
line  in  Figure  4  compared  to  those  derived  from  the  tele¬ 
graphist’s  equations.  The  metals  form  ohmic  contacts  with 
the  p-type  silicon  substrate.  The  simulated  data  labeled 
“full  wave”  were  computed  with  full  wave  electromagnetics, 
and  those  labeled  “quasi”  used  the  quasi-static  field  approx¬ 
imation. 


4.  RESULTS  AND  DISCUSSION 

To  test  the  simulation,  the  simple  lossy  transmission 
line  shown  in  Figure  4  is  considered.  It  consists  of  two 


Fig.  4:  A  parallel  conductor  transmission  line  with  a  doped 
semiconductor  serving  as  a  lossy  substrate  and  a  50  ter¬ 
mination  at  the  output. 

ideal  metal  contacts  on  a  doped  semiconductor  substrate 
terminated  with  a  50  impedance.  For  the  first  test  case, 
the  metals  form  ohmic  contacts  with  the  substrate,  which 


wave  data  compare  favorably  with  the  analytic  solutions, 
with  the  increasing  error  at  higher  frequencies  most  likely 
due  to  decreasing  Fourier  components  of  the  input  signal. 
A  flatter  input  spectrum  decreases  the  error.  Figure  4  also 
shows  very  poor  agreement  between  the  quasi-static  and 
analytic  solutions.  This  is  interesting  because  full  wave  ef¬ 
fects  are  typically  thought  of  in  terms  of  interference  that 
occurs  when  feature  sizes  are  comparable  to  wavelength. 
When  feature  sizes  are  much  smaller  than  the  wavelength, 
the  quasi-static  field  approximation  is  typically  considered 
appropriate.  However,  the  transmission  line  in  Figure  4  is 
an  order  of  magnitude  shorter  than  the  smallest  signal  wave¬ 
length  considered,  and  yet  the  quasi-static  approximation  is 
highly  inaccurate. 

The  electrically  short  transmission  line  exhibits  the  sig¬ 
nificant  full  wave  effect  shown  in  Figure  5  because  quasi¬ 
static  fields  do  not  store  energy.  When  the  Gaussian  input 
signal  is  applied,  energy  is  injected  into  the  line,  and  that 
energy  will  persist  for  as  long  as  it  takes  the  terminating 
impedance  and  the  mobile  charges  in  the  substrate  to  dis¬ 
sipate  it.  The  energy  is  stored  in  the  form  of  rotational 
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electric  and  magnetic  fields,  determined  by  Ampere’s  and 
Faraday’s  laws.  Since  the  quasi-static  approximation  does 


Fig.  6:  Simulated  output  currents  for  the  lossy  transmission 
line  in  Figure  4  excited  by  a  Gaussian  input  voltage  signal. 
The  currents  labeled  “full  wave”  were  computed  with  full 
wave  electromagnetics,  and  those  labeled  “quasi”  used  the 
quasi-static  field  approximation. 

not  consider  these  fields,  the  current  delivered  to  the  output 
impedance  simply  follows  the  input  voltage  and  is  unable 
to  represent  the  currents  that  persist  after  the  input  signal 
goes  to  zero,  as  shown  in  Figure  6. 

A  nonlinear  example,  in  which  the  top  metal  in  Fig¬ 
ure  4  forms  a  Schottky  contact  with  a  doped  (1015  cm~3  n- 
type)  GaAs  substrate,  was  also  considered.  The  rectifying 
Schottky  contact  blocks  current  for  negative  applied  volt¬ 
ages  but  allows  the  substrate  to  conduct  significant  currents 
under  sufficiently  large  positive  voltages.  It  also  produces 
a  very  narrow  depletion  region  beneath  the  top  metal  that 
greatly  increases  the  transmission  line’s  capacitance.  Fig¬ 
ure  7  shows  S-parameter  computed  around  quiescent  bias 
points  of  0  V  (nonconducting)  and  2  V  (highly  conducting). 
The  nonlinear  properties  of  this  transmission  line  make  its 
frequency  characteristics  strong  functions  of  its  bias  condi¬ 
tion,  and  like  the  previous  example,  high  frequency  solutions 
require  full  wave  electromagnetics. 

4.  CONCLUSIONS 

Maxwell’s  field  equations  and  the  thermodynamics  of 
ideal  Fermi  gases  of  charged  particles  form  a  complete  and 
physically  self-consistent  classical  model  of  the  electromag¬ 
netic  interaction.  The  remarkable  compatibility  of  the  vec¬ 
tor  calculus  of  the  fields  and  the  statistical  mechanics  of  the 
charges  is  not  coincidental.  Since  Maxwell  developed  the 
field  theory  while  he  was  also  working  on  the  kinetic  theory 
of  gases,  it  is  likely  the  two  theories  were  intended  to  work 
together.  Although  the  quantum  mechanics  of  fermions  re¬ 
quired  the  Maxwell-Boltzmann  distribution  of  the  charged 
particles  to  be  replaced  with  the  more  general  Fermi-Dirac 
distribution,  the  essential  physics  is  the  same.  The  field 
theory  has  not  required  any  modification  from  its  original 
form. 


Fig.  7:  Simulated  S-parameters  for  the  lossy  transmission 
line  in  Figure  4  with  the  top  metal  forming  a  Schottky  con¬ 
tact  with  an  n-type  GaAs  substrate.  The  voltage  labels 
show  the  quiescent  biases.  Curves  labeled  “full  wave”  were 
computed  with  full  wave  electromagnetics,  and  those  la¬ 
beled  “quasi”  used  the  quasi-static  field  approximation. 


Although  the  necessary  pieces  of  a  complete  classical 
theory  of  the  electromagnetic  interaction  have  been  in  place 
for  a  long  time,  the  model  had  not  been  computed  self- 
consistently.  The  most  commonly  developed  techniques 
for  treating  the  vector  field  components  of  the  model  have 
not  been  well  suited  for  highly  nonlinear  charge  transport, 
and  the  techniques  developed  for  nonlinear  charge  trans¬ 
port  have  not  been  generalized  to  treat  vector  fields  self- 
consistently.  This  paper  has  presented  a  new  discretization 
scheme  called  DVSI  that  overcomes  these  obstacles  by  rep¬ 
resenting  both  the  curls  of  vector  fields  and  the  divergences 
of  fluxes  in  a  self-consistent  manner. 

DVSI  permits  the  fully  coupled  solution  of  the  classi¬ 
cal  electromagnetic  interaction  for  highly  nonlinear  systems 
from  dc  to  arbitrarily  high  frequencies.  It  was  used  to  sim¬ 
ulate  lossy  transmission  lines  that  exhibit  different  degrees 
of  nonlinearity.  A  highly  linear  case  showed  the  simulated 
results  to  agree  well  with  analytic  solutions.  This  simple  re¬ 
sult  also  demonstrated  that  full  wave  electromagnetics  can 
be  important  even  when  devices  features  are  much  smaller 
than  signal  wavelengths.  This  is  true  when  the  device  is 
able  to  store  electromagnetic  energy.  The  energy  is  stored 
in  the  form  of  rotational  fields  determined  by  Ampere’s  and 
Faraday’s  laws.  Additional  simulations  of  nonlinear  trans¬ 
mission  lines  showed  that  the  high  frequency  characteristics 
strongly  depend  on  the  dc  quiescent  bias  condition.  Further 
potential  uses  for  the  code  include  the  full  wave  simulation 
of  high  speed  microelectronics  and  optoelectronics,  which 
will  be  the  subjects  of  future  work. 
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